You can download these files by clicking on the links
Working with Luminosity
Note that these are large files: Every zip file will range from 355 MB to 411.1 MB
The files are cloud-free composites
In cases where two satellites were collecting data - two composites were produced.
Working with Luminosity
Working with Luminosity
Note that these are large files: Every zip file will range from 355 MB to 411.1 MB
The files are cloud-free composites
In cases where two satellites were collecting data - two composites were produced.
The products are 30 arc second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude.
Each composite set is named with the satellite and the year (F121995 is from DMSP satellite number F12 for the year 1995)
Working with Luminosity
Three types of images are available:
F1?YYYY_v4b_cf_cvg.tif: Cloud-free coverages tally the total number of observations that went into each 30 arc second grid cell
F1?YYYY_v4b_avg_vis.tif Raw avg_vis contains the average of the visible band digital number values with no further filtering.
F1?YYYY_v4b_stable_lights.avg_vis.tif The cleaned up avg_vis contains the lights from cities, towns, and other sites with persistent lighting, including gas flares.
Working with Luminosity
The DMSP/OLS resolution is approximately 1000m
It comprises six different DMSP satellites:
F10 (1992–1994)
F12 (1994–1999)
F14 (1997–2003)
F15 (2000–2007)
F16 (2004–2009)
F18 (2010–2013)
Working with Luminosity
After downloading, your folder should look like below:
Working with Luminosity
You download all the relevant files for this exercise below:
These are the steps to have the data ready for plotting in ggplot
#Step1: Read country bordersborders <-st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet =TRUE)#Step2: Creating a folder calledtarget_dir <-"./data/luminosity/F101992"#Step3: Create the target directory if it doesn't existif (!dir.exists(target_dir)) {dir.create(target_dir)}#Step4: Unzip the fileuntar("./data/luminosity/F101992.v4.tar", exdir = target_dir)#Step5: Unzip the file furtherR.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove =FALSE, overwrite=TRUE)#Step6: Read the rasterf1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")#Step7: Make sure the two files have the same CRSst_crs(f1992)<-st_crs(borders)#Step8: Reduce the resolution of the rasterf1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx =0.7))#Step9: Rename the rasternames(f1992b)<-"lights_1992"fig<-ggplot()+geom_stars(data=f1992b)+scale_fill_gradientn(colours=c("black","white"))+#scale_fill_viridis_c(option = "magma",begin = 0.1)+geom_sf(data=borders, linewidth =0.1, fill =NA, color ="white", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())fig
Luminosity in 1992
Luminosity in 1992 in Europe
Raster Values and Properties
The print method gives a summary of the raster’s properties.
print(f1992b)
stars object with 2 dimensions and 1 attribute
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
lights_1992 0 2 3 3.384738 4 63 515
dimension(s):
from to offset delta refsys x/y
x 1 515 -180 0.7 WGS 84 [x]
y 1 201 75 -0.7 WGS 84 [y]
The class function allows us to see the type of object that we’re dealing with
class(f1992b)
[1] "stars"
Raster Values and Properties
Fundamentally, a stars object is a collection of matrix or array objects.
This is along with additional properties of the dimensions, such as dimension names, Coordinate Reference Systems (CRS).
We can display the structure of a stars object using str.
str(f1992b)
List of 1
$ lights_1992: num [1:515, 1:201] 0 0 0 0 0 0 0 0 0 0 ...
- attr(*, "dimensions")=List of 2
..$ x:List of 7
.. ..$ from : num 1
.. ..$ to : num 515
.. ..$ offset: num -180
.. ..$ delta : num 0.7
.. ..$ refsys:List of 2
.. .. ..$ input: chr "WGS 84"
.. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
.. .. ..- attr(*, "class")= chr "crs"
.. ..$ point : logi NA
.. ..$ values: NULL
.. ..- attr(*, "class")= chr "dimension"
..$ y:List of 7
.. ..$ from : num 1
.. ..$ to : num 201
.. ..$ offset: num 75
.. ..$ delta : num -0.7
.. ..$ refsys:List of 2
.. .. ..$ input: chr "WGS 84"
.. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
.. .. ..- attr(*, "class")= chr "crs"
.. ..$ point : logi NA
.. ..$ values: NULL
.. ..- attr(*, "class")= chr "dimension"
..- attr(*, "raster")=List of 4
.. ..$ affine : num [1:2] 0 0
.. ..$ dimensions : chr [1:2] "x" "y"
.. ..$ curvilinear: logi FALSE
.. ..$ blocksizes : NULL
.. ..- attr(*, "class")= chr "stars_raster"
..- attr(*, "class")= chr "dimensions"
- attr(*, "class")= chr "stars"
Raster Values and Properties
From this, we see that we have the f1992b object of length 1.
The 2nd line in the str printout specifies the name and contents of the first (and only) item in the list, namely its type, dimensions, and a sample of the first few values.
In our case, the first item is lights_1992, and it is a matrix with 515 rows and 201 columns
Raster Attributes and Values
To access some of the attributes of our objects, we can type:
names(f1992b)
[1] "lights_1992"
We can change names easily in the following wat:
names(f1992b)<-"lum_1992"names(f1992b)
[1] "lum_1992"
Accessing Attributes by Name
We can access our raster’s attributes by using the list access methods
A histogram will help us identify the raster values distribution
#Turning the object into a dataframedf<-st_as_sf(f1992b[1], as_points =FALSE, merge =FALSE)#Creating a histogramggplot(df, aes(x = lum_1992)) +geom_histogram(binwidth =10, color ="white") +ggtitle("Histogram of lum_1992") +xlab("lum_1992") +ylab("Frequency")+theme_bw()
Accessing Raster Values
We thus see that most of values around the world are 0.
This is because 71 percent of the Earth’s surface is water-covered.
The important part is that by transforming the stars object into a dataframe (grid), we can perform the usual operations dedicated to dataframes.
#Step1: Reading the central polygonrelev_poly<-st_read(dsn="./data/rome_center/central_polygon.shp", quiet = T)#Step2: Reading roadsgis_osm_roads <-st_read(dsn="./data/selection.gdb", layer="gis_osm_roads", quiet = T)#Step3: Extracting the bix from the central polygonext2<-st_bbox(relev_poly)#Step4: Ensuring the same crs for the polygon as for the rasterst_crs(relev_poly)<-st_crs(f1992)#Step5: Cropping the rasterrc <-st_crop(f1992, relev_poly)#Step6: Mapping everythingggplot()+geom_stars(data = rc)+geom_sf(data = relev_poly, color="white", fill=NA)+geom_sf(data = gis_osm_roads, color="white", linewidth=0.05)+theme_bw()+scale_fill_gradientn(colours=c("black","white"))
We can easily replace the values inside the raster (grid) by using simple dataframe operations.
For example, let’s replace all the values that have 63 with NA:
Let us now make the NA grids transparent
Writing raster to file
We can of course write the grid as a shapefile for we can turn it back to a raster.
rc3_raster =st_rasterize(rc3)class(rc3_raster)
[1] "stars"
class(rc3)
[1] "sf" "data.frame"
#Step1: Creating a folder calledtarget_dir <-"./data/output"#Step2: Create the target directory if it doesn't existif (!dir.exists(target_dir)) {dir.create(target_dir)}#Step3: Writing the rasterwrite_stars(rc3_raster, './data/output/rc3_raster.tif', update =TRUE)#Step4: Writing the shape filesst_write(rc3, "./data/output/rc3_grid.shp", delete_dsn =TRUE)
Deleting source `./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing layer `rc3_grid' to data source
`./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing 535 features with 1 fields and geometry type Polygon.
Working with More Complex Rasters
Working with lights was simple
In many other contexts, we can have more complex rasters
Temperature data from CEDA is an example of a raster that contains time series